Python notebook¶

Met deze Notebook genereren we de volgende output:

• vertical_resistance_layer1: een nieuwe verticale weerstand voor de deklaag
• winter\peil_laag#_#.idf: nieuwe peilen voor de winter in de verschillende (sub) lagen
• zomer\peil_laag#_#.idf: nieuwe peilen voor de zomer in de verschillende (sub) lagen

Een aantal relevante tussenproducten:

• peilen.gpkg: de geaggregeerde vaste/zomer/winter en flexibele peilen naar 1 zomer en 1 winterpeil, wanneer deze kan worden berekend
• winterpeil_ruw.tif, winterpeil.tif, winterpeil.idf: het ruwe en geinterpoleerde resultaat van de verrastering van de kolom winterpeil in peilen.gpkg
• zomerpeil_ruw.tif, zomerpeil.tif, zomerpeil.idf: het ruwe en geinterpoleerde resultaat van de verrastering van de kolom zomerpeil in peilen.gpkg

In deze Notebook worden alle tussen en eindproducten stap-voor-stap gegenereerd met duiding en Python-code.

In [1]:

try: # to bypass ImportError at rasterio.__init__ from rasterio.__version ....
    import rasterio
except ImportError:
    from osgeo import gdal
    import rasterio

from dask import array # If you don't do this imod.idf.open() (version 0.9.0) will give an AttributeError: module 'dask' has no attribute 'array'
import geopandas as gpd
import imod
import numpy as np
from pathlib import Path
import shutil
from rasterio.features import rasterize
from rasterio.fill import fillnodata
import xarray as xr

from nhi_data import Ondergrond

import warnings
warnings.filterwarnings("ignore") # To avoid logging RuntimeWarning in rasterio.fill.fillnodata

try: # to bypass ImportError at rasterio.__init__ from rasterio.__version .... import rasterio except ImportError: from osgeo import gdal import rasterio from dask import array # If you don't do this imod.idf.open() (version 0.9.0) will give an AttributeError: module 'dask' has no attribute 'array' import geopandas as gpd import imod import numpy as np from pathlib import Path import shutil from rasterio.features import rasterize from rasterio.fill import fillnodata import xarray as xr from nhi_data import Ondergrond import warnings warnings.filterwarnings("ignore") # To avoid logging RuntimeWarning in rasterio.fill.fillnodata

Specificaties iMODFLOW inlezen¶

We lezen de specificatie van het iMODFLOW model in. We gebruiken daarvoor de laag PEIL_LAAG1_1.IDF, waarvan we de ruimtelijke attributen overnemen.

In [2]:

oppervlaktewater_dir = Path(r"d:\projecten\D2304.HDSR_TEO_debieten\02.model\data\hydromedah\oppervlaktewater")

template_idf = oppervlaktewater_dir.joinpath("winter", "PEIL_LAAG1_1.IDF")
template_da = imod.idf.open(template_idf)

shape = template_da.shape
height, width = template_da.shape
transform = imod.util.transform(template_da)

_, xmin, xmax, _, ymin, ymax = imod.util.spatial_reference(template_da)

oppervlaktewater_dir = Path(r"d:\projecten\D2304.HDSR_TEO_debieten\02.model\data\hydromedah\oppervlaktewater") template_idf = oppervlaktewater_dir.joinpath("winter", "PEIL_LAAG1_1.IDF") template_da = imod.idf.open(template_idf) shape = template_da.shape height, width = template_da.shape transform = imod.util.transform(template_da) _, xmin, xmax, _, ymin, ymax = imod.util.spatial_reference(template_da)

Vervangen deklaagweerstand¶

Hier vervangen we de verticale weerstand van de deklaag uit het LHM4.1

Inladen NHI ondergrond¶

We maken een Ondergrond-object, voor de bounding-box van het HDSR MODFLOW model. We zoeken de juiste laag op de NHI-geoserver

In [3]:

ondergrond = Ondergrond(bbox=(xmin,ymin,xmax,ymax))
layer = ondergrond.get_layers(filter="resistance_layer1")[0]
print(layer)

ondergrond = Ondergrond(bbox=(xmin,ymin,xmax,ymax)) layer = ondergrond.get_layers(filter="resistance_layer1")[0] print(layer)

ondergrond2:vertical_resistance_layer1_lhm41

Download en conversie naar GeoTIFF en IDF¶

Met dit object kunnen we de laag downloaden en opslaan als GeoTiff en IDF

In [4]:

ondergrond.download_layer(layer, "vertical_resistance_layer1.tif")
ondergrond.download_layer(layer, "vertical_resistance_layer1.idf")

ondergrond.download_layer(layer, "vertical_resistance_layer1.tif") ondergrond.download_layer(layer, "vertical_resistance_layer1.idf")

Vervangen peil en peilgebieden¶

Hier vervangen we de oppervlaktewaterpeilen van het iMODFLOW model

Inladen peilgebieden¶

De peilgebieden zijn op 04-07-2023 gedownload (gegenereerd op 23 jun 2023 11:04) via de HDSR Infovijver als GeoJSON.

In [5]:

peilengebieden_gdf = gpd.read_file("Peilgebieden.geojson", engine="pyogrio")
peilengebieden_gdf.explore()

peilengebieden_gdf = gpd.read_file("Peilgebieden.geojson", engine="pyogrio") peilengebieden_gdf.explore()

Out[5]:

Make this Notebook Trusted to load map: File -> Trust Notebook

Aggregeren peilen naar zomer- en winterpeil¶

We lezen de geometrieen in van de peilgebieden en kopieren het attribuut SOORTSTREEFPEIL (soort_streefpeil) met de volgende enumeratie:

• 1: vast peil
• 2: flexibel peil
• 3: zomer- en winterpeil
• 4, 5, 99: geen peil vastgesteld/peil onbekend

We willen de volgende kolommen vullen:

• code: waterschapscode
• nodata: geen zomer- winterpeil bekend = True, anders False
• betrouwbaar_peil: winterpeil <= zomerpeil = True, anders False
• zomerpeil
• winterpeil

We projecteren de geometriën naar Rijksdriehoekstelsel (epsg 28992) en vullen de kolommen soort_streefpeil, code uit de oorspronkelijke data en zetten de waarden van de kolom nodata op False.

In [6]:

peilen_gdf = gpd.GeoDataFrame(peilengebieden_gdf["geometry"])
peilen_gdf.to_crs(28992, inplace=True)
peilen_gdf["soort_streefpeil"] = peilengebieden_gdf["SOORTSTREEFPEIL"]
peilen_gdf["code"] = peilengebieden_gdf["CODE"]
peilen_gdf["nodata"] = False

peilen_gdf = gpd.GeoDataFrame(peilengebieden_gdf["geometry"]) peilen_gdf.to_crs(28992, inplace=True) peilen_gdf["soort_streefpeil"] = peilengebieden_gdf["SOORTSTREEFPEIL"] peilen_gdf["code"] = peilengebieden_gdf["CODE"] peilen_gdf["nodata"] = False

Vast Peil¶

Wanneer VASTPEIL == -999 of niet ingevuld, dan geven we dit aan in nodata. Voor de peilgebieden waarvan gegevens bekend zijn geldt zomerpeil = winterpeil = vastpeil

In [7]:

peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] == 1

# VASTPEIL = -999 komt voor, dit maskeren we en leggen we vast als no data
data_mask = (peilengebieden_gdf["VASTPEIL"] != -999) & ~peilengebieden_gdf["VASTPEIL"].isna()
peilen_gdf.loc[(peil_mask &  ~data_mask), ["nodata"]] = True

# vullen zomerpeil en winterpeil
peilen_gdf.loc[(peil_mask & data_mask), ["zomerpeil"]] = peilengebieden_gdf[(peil_mask & data_mask)]["VASTPEIL"]
peilen_gdf.loc[(peil_mask & data_mask), ["winterpeil"]] = peilengebieden_gdf[(peil_mask & data_mask)]["VASTPEIL"]

peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] == 1 # VASTPEIL = -999 komt voor, dit maskeren we en leggen we vast als no data data_mask = (peilengebieden_gdf["VASTPEIL"] != -999) & ~peilengebieden_gdf["VASTPEIL"].isna() peilen_gdf.loc[(peil_mask & ~data_mask), ["nodata"]] = True # vullen zomerpeil en winterpeil peilen_gdf.loc[(peil_mask & data_mask), ["zomerpeil"]] = peilengebieden_gdf[(peil_mask & data_mask)]["VASTPEIL"] peilen_gdf.loc[(peil_mask & data_mask), ["winterpeil"]] = peilengebieden_gdf[(peil_mask & data_mask)]["VASTPEIL"]

Flexibel peil¶

In zowel ONDERPEIL als BOVENPEIL komt -999 voor, maar niet gelijktijdig. In deze gevallen stellen we ONDERPEIL gelijk aan BOVENPEIL of vise-versa. Vervolgens brekenen we zomerpeil en winterpeil:

zomerpeil = winterpeil = (ONDERPEIL + BOVENPEIL) / 2

In [8]:

peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] == 2

# bij -999 in bovenpeil stellen we bovenpeil gelijk aan onderpeil
nodata_mask = peil_mask & (peilengebieden_gdf["BOVENPEIL"] == -999)
peilengebieden_gdf.loc[nodata_mask, ["BOVENPEIL"]] = peilengebieden_gdf[nodata_mask]["ONDERPEIL"]

# bij -999 in onderpeil stellen we onderpeil gelijk aan bovenpeil
nodata_mask = peil_mask & (peilengebieden_gdf["ONDERPEIL"] == -999)
peilengebieden_gdf.loc[nodata_mask, ["ONDERPEIL"]] = peilengebieden_gdf[nodata_mask]["BOVENPEIL"]

peilen_gdf.loc[peil_mask, ["zomerpeil"]] = (peilengebieden_gdf[peil_mask]["ONDERPEIL"] + peilengebieden_gdf[peil_mask]["BOVENPEIL"]) / 2
peilen_gdf.loc[peil_mask, ["winterpeil"]] = peilen_gdf[peil_mask]["zomerpeil"]

peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] == 2 # bij -999 in bovenpeil stellen we bovenpeil gelijk aan onderpeil nodata_mask = peil_mask & (peilengebieden_gdf["BOVENPEIL"] == -999) peilengebieden_gdf.loc[nodata_mask, ["BOVENPEIL"]] = peilengebieden_gdf[nodata_mask]["ONDERPEIL"] # bij -999 in onderpeil stellen we onderpeil gelijk aan bovenpeil nodata_mask = peil_mask & (peilengebieden_gdf["ONDERPEIL"] == -999) peilengebieden_gdf.loc[nodata_mask, ["ONDERPEIL"]] = peilengebieden_gdf[nodata_mask]["BOVENPEIL"] peilen_gdf.loc[peil_mask, ["zomerpeil"]] = (peilengebieden_gdf[peil_mask]["ONDERPEIL"] + peilengebieden_gdf[peil_mask]["BOVENPEIL"]) / 2 peilen_gdf.loc[peil_mask, ["winterpeil"]] = peilen_gdf[peil_mask]["zomerpeil"]

Zomer- en winterpeil¶

Hier komen geen ontbrekende, danwel -999 waarden voor, dus deze nemen we direct over

In [9]:

peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] == 3

peilen_gdf.loc[peil_mask, ["zomerpeil"]] = peilengebieden_gdf[peil_mask]["ZOMERPEIL"]
peilen_gdf.loc[peil_mask, ["winterpeil"]] = peilengebieden_gdf[peil_mask]["WINTERPEIL"]

peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] == 3 peilen_gdf.loc[peil_mask, ["zomerpeil"]] = peilengebieden_gdf[peil_mask]["ZOMERPEIL"] peilen_gdf.loc[peil_mask, ["winterpeil"]] = peilengebieden_gdf[peil_mask]["WINTERPEIL"]

Onbrekende peilen¶

Wanneer peilen ontbreken, geven we dat aan door nodata op True te zetten

In [10]:

peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] >= 4
peilen_gdf.loc[peil_mask, ["nodata"]] = True

peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] >= 4 peilen_gdf.loc[peil_mask, ["nodata"]] = True

Indicatie betrouwbaarheid¶

Voor alle gebieden met zomerpeil en winterpeil controleren we of winterpeil kleiner of gelijk aan zomerpeil is.

In [11]:

data_mask = ~peilen_gdf["nodata"]
peilen_gdf.loc[data_mask, ["peil_betrouwbaar"]] =  peilen_gdf[data_mask]["winterpeil"] <= peilen_gdf[data_mask]["zomerpeil"]

data_mask = ~peilen_gdf["nodata"] peilen_gdf.loc[data_mask, ["peil_betrouwbaar"]] = peilen_gdf[data_mask]["winterpeil"] <= peilen_gdf[data_mask]["zomerpeil"]

Wegschrijven¶

In [12]:

peilen_gdf.to_file("peilen.gpkg")

peilen_gdf.to_file("peilen.gpkg")

Conversie naar imod¶

Het resultaat in peilen_gdf wordt weggeschreven naar idf-bestanden. Van de verrasterde peilen schrijven we tevens GeoTifs weg, zodat deze in GIS kunnen worden geladen.

In [13]:

crs = peilen_gdf.crs

profile={ # profiel voor GeoTifs
    "driver": "GTiff",
    "dtype":rasterio.dtypes.float32,
    "nodata": -999,
    "width" :width,
    "height": height,
    "count": 1,
    "crs": crs,
    "transform": transform}

crs = peilen_gdf.crs profile={ # profiel voor GeoTifs "driver": "GTiff", "dtype":rasterio.dtypes.float32, "nodata": -999, "width" :width, "height": height, "count": 1, "crs": crs, "transform": transform}

Data mask¶

We zorgen dat de Utrechtse Heuvelrug (PG2179) en Nederrijn en Lek (PG2112) uit de peilen-rasters blijven. We houden daar de oorspronkelijke waarden. Het deel waar we peilen voor verrasteren schrijven we weg als mask.tif.

In [14]:

shapes = (
    (geom, 1)
    for geom in peilen_gdf[~peilen_gdf["code"].isin(["PG2179", "PG2112"])]["geometry"]
    )

mask_array = rasterize(
    shapes=shapes,
    out_shape=shape,
    transform=transform,
    all_touched=True)

with rasterio.open("mask.tif", "w+", **profile) as dst:
    dst.write(mask_array, 1)

shapes = ( (geom, 1) for geom in peilen_gdf[~peilen_gdf["code"].isin(["PG2179", "PG2112"])]["geometry"] ) mask_array = rasterize( shapes=shapes, out_shape=shape, transform=transform, all_touched=True) with rasterio.open("mask.tif", "w+", **profile) as dst: dst.write(mask_array, 1)

Verwerking in IDFs¶

Voor de seizoenen zomer en winter voeren we de volgende stappen uit:

1. Verrasteren van de zomer en winterpeilen waar deze data bevatten. We schrijven deze weg in zomerpeil_ruw.tif en winterpeil_ruw.tif
2. Het opvullen van nodata cellen binnen het mask door middel van interpolatie. We schrijven deze weg als zomerpeil.tif en winterpeil.tif
3. Het overschrijven van de waarden in alle PEIL_LAAG*.IDF bestanden uit het oorsprokelijke model met de regels:
   • nodata in de oorspronkelijke peilen blijft nodatain de nieuwe peilen
   • Wanneer er een er een bodemhoogte is wordt het nieuwe peil gelijk of hoger dan de bodemhoogte

In [15]:

for seizoen in ["zomer", "winter"]:
    # 1. Verrasteren van de zomer en winterpeilen waar deze data bevatten. We schrijven deze weg in `zomerpeil_origineel.tif` en `winterpeil_origineel.tif`
    shapes = (
        (geom, value)
        for geom, value in zip(peilen_gdf[~peilen_gdf.nodata]["geometry"], peilen_gdf[~peilen_gdf.nodata][f"{seizoen}peil"])
        )
    
    peil_array = rasterize(
        shapes=shapes,
        out_shape=shape,
        fill=-999,
        transform=transform,
        all_touched=False)
    
    with rasterio.open(f"{seizoen}peil_ruw.tif", "w+", **profile) as dst:
        dst.write(peil_array, 1)

    # 2. Het opvullen van `nodata` cellen binnen het mask door middel van interpolatie. We schrijven deze weg als `zomerpeil.tif` en `winterpeil.tif
    mask = np.where((peil_array == -999) & (mask_array == 1), 0, 1)
    mask = np.where(mask_array == 0, np.NaN, mask)
    peil_array = fillnodata(
        peil_array,
        mask=mask,
        max_search_distance=20,
        smoothing_iterations=0
        )
    
    peil_array = np.where(mask_array == 1, peil_array, -999)
    with rasterio.open(f"{seizoen}peil.tif", "w+", **profile) as dst:
        dst.write(peil_array, 1)
    
    dataset = imod.rasterio.open(f"{seizoen}peil.tif")
    imod.idf.save(f"{seizoen}peil", dataset)

    # 3. Het overschrijven van de waarden in alle `PEIL_LAAG*.IDF` bestanden
    out_dir = Path(seizoen)
    if out_dir.exists():
        shutil.rmtree(out_dir)
    out_dir.mkdir()
    datasets = imod.idf.open_dataset(oppervlaktewater_dir.joinpath(seizoen,"PEIL_LAAG*.IDF"))
    
    for laag, peil_da in datasets.items():

        # - `nodata` in de oorspronkelijke peilen blijft `nodata`in de nieuwe peilen
        peil_da = xr.where(
            np.isnan(peil_da),
            peil_da,
            np.where(peil_array == -999, peil_da, peil_array)
            )
        
        # - Wanneer er een er een bodemhoogte is wordt het nieuwe peil gelijk of hoger dan de bodemhoogte
        bodemhoogte_name = f"BODEMHOOGTE_LAAG{laag[len(laag)-3:]}"
        bodemhoogte_da = imod.idf.open(
            oppervlaktewater_dir.joinpath(seizoen,f"{bodemhoogte_name}.IDF")
            )
        peil_da = xr.where(
            np.isnan(peil_da),
            peil_da,
            xr.where(
                np.isnan(bodemhoogte_da),
                peil_da,
                xr.where(
                    bodemhoogte_da < peil_da,
                    peil_da,
                    bodemhoogte_da
                    ) # zowel bodemhoogte als peil zijn niet NaN
                )
            )
        # Wegschrijven IDF
        imod.idf.write(out_dir / f"{laag}.idf", peil_da)

for seizoen in ["zomer", "winter"]: # 1. Verrasteren van de zomer en winterpeilen waar deze data bevatten. We schrijven deze weg in `zomerpeil_origineel.tif` en `winterpeil_origineel.tif` shapes = ( (geom, value) for geom, value in zip(peilen_gdf[~peilen_gdf.nodata]["geometry"], peilen_gdf[~peilen_gdf.nodata][f"{seizoen}peil"]) ) peil_array = rasterize( shapes=shapes, out_shape=shape, fill=-999, transform=transform, all_touched=False) with rasterio.open(f"{seizoen}peil_ruw.tif", "w+", **profile) as dst: dst.write(peil_array, 1) # 2. Het opvullen van `nodata` cellen binnen het mask door middel van interpolatie. We schrijven deze weg als `zomerpeil.tif` en `winterpeil.tif mask = np.where((peil_array == -999) & (mask_array == 1), 0, 1) mask = np.where(mask_array == 0, np.NaN, mask) peil_array = fillnodata( peil_array, mask=mask, max_search_distance=20, smoothing_iterations=0 ) peil_array = np.where(mask_array == 1, peil_array, -999) with rasterio.open(f"{seizoen}peil.tif", "w+", **profile) as dst: dst.write(peil_array, 1) dataset = imod.rasterio.open(f"{seizoen}peil.tif") imod.idf.save(f"{seizoen}peil", dataset) # 3. Het overschrijven van de waarden in alle `PEIL_LAAG*.IDF` bestanden out_dir = Path(seizoen) if out_dir.exists(): shutil.rmtree(out_dir) out_dir.mkdir() datasets = imod.idf.open_dataset(oppervlaktewater_dir.joinpath(seizoen,"PEIL_LAAG*.IDF")) for laag, peil_da in datasets.items(): # - `nodata` in de oorspronkelijke peilen blijft `nodata`in de nieuwe peilen peil_da = xr.where( np.isnan(peil_da), peil_da, np.where(peil_array == -999, peil_da, peil_array) ) # - Wanneer er een er een bodemhoogte is wordt het nieuwe peil gelijk of hoger dan de bodemhoogte bodemhoogte_name = f"BODEMHOOGTE_LAAG{laag[len(laag)-3:]}" bodemhoogte_da = imod.idf.open( oppervlaktewater_dir.joinpath(seizoen,f"{bodemhoogte_name}.IDF") ) peil_da = xr.where( np.isnan(peil_da), peil_da, xr.where( np.isnan(bodemhoogte_da), peil_da, xr.where( bodemhoogte_da < peil_da, peil_da, bodemhoogte_da ) # zowel bodemhoogte als peil zijn niet NaN ) ) # Wegschrijven IDF imod.idf.write(out_dir / f"{laag}.idf", peil_da)